#ifndef SPHERICITY_H
#define SPHERICITY_H

#include "TMath.h"
#include "TMatrix.h"
#include "TRandom.h"
#include "TVector3.h"
#include "TLorentzVector.h"
#include "TObjArray.h"
#include "TClass.h"

/*
class TMatrix;
class TVector3;
class TLorentzVector;
*/
class Sphericity: public TObject{
	

public:
// Constructor.
	Sphericity() {

// Minimum number of particles to perform study.
NSTUDYMIN     = 2;

// Maximum number of times that an error warning will be printed.
TIMESTOPRINT  = 1;

// Assign mimimum squared momentum in weight to avoid division by zero. 
P2MIN         = 1e-20;

// Second eigenvalue not too low or not possible to find eigenvectors.
EIGENVALUEMIN = 1e-10;

}
	~Sphericity() {}
void setParams (double powerIn = 2.) {
	power = powerIn;
	nFew=0;
	 nBack=0;
     powerInt = 0; 
    if (fabs(power - 1.) < 0.01) powerInt = 1;
    if (fabs(power - 2.) < 0.01) powerInt = 2; 
    powerMod = 0.5 * power - 1.;

return;
}
  
  // Analyze event.
  bool analyze(TObjArray* e, double powerIn = 2.)
{


setParams(powerIn );

  // Initial values, tensor and counters zero.
  eVal1 = eVal2 = eVal3 = 0.;
  eVec1 = eVec2 = eVec3 = 0.;
  double tt[4][4];
  for (int j = 1; j < 4; ++j) 
  for (int k = j; k < 4; ++k) tt[j][k] = 0.;
  int nStudy = 0;
  double denom = 0.;

  // Loop over desired particles in the event.
  for (int i = 0; i < e->GetEntries(); ++i) {
 TObject* o = e->At(i);
 TString nam(o->IsA()->GetName());
  double pNow[4];

 if (nam.Contains("TVector3")) {
    pNow[1] = ((TVector3 *) o)->X();
    pNow[2] = ((TVector3 *) o)->Y();
    pNow[3] = ((TVector3 *) o)->Z();
} else if (nam.Contains("TLorentzVector")) {
    pNow[1] = ((TLorentzVector *) o)->X();
    pNow[2] = ((TLorentzVector *) o)->Y();
    pNow[3] = ((TLorentzVector *) o)->Z();

	}
  

    ++nStudy;

    // Calculate matrix to be diagonalized. Special cases for speed.
    
    double p2Now = pNow[1]*pNow[1] + pNow[2]*pNow[2] + pNow[3]*pNow[3];
    double pWeight = 1.;
    if (powerInt == 1) pWeight = 1. / sqrt(max(P2MIN, p2Now));
    else if (powerInt == 0) pWeight = pow( max(P2MIN, p2Now), powerMod);
    for (int j = 1; j < 4; ++j)   
    for (int k = j; k < 4; ++k) tt[j][k] += pWeight * pNow[j] * pNow[k];
    denom += pWeight * p2Now;
  
}
  // Very low multiplicities (0 or 1) not considered.
  if (nStudy < NSTUDYMIN) {
    if (nFew < TIMESTOPRINT) 
    ++nFew;
    return false;
  }

  // Normalize tensor to trace = 1.
  for (int j = 1; j < 4; ++j) 
  for (int k = j; k < 4; ++k) tt[j][k] /= denom;
 
  // Find eigenvalues to matrix (third degree equation).
  double qCoef = ( tt[1][1] * tt[2][2] + tt[1][1] * tt[3][3] 
    + tt[2][2] * tt[3][3] - pow2(tt[1][2]) - pow2(tt[1][3]) 
    - pow2(tt[2][3]) ) / 3. - 1./9.;
  double qCoefRt = sqrt( -qCoef);
  double rCoef = -0.5 * ( qCoef + 1./9. + tt[1][1] * pow2(tt[2][3]) 
    + tt[2][2] * pow2(tt[1][3]) + tt[3][3] * pow2(tt[1][2]) 
    - tt[1][1] * tt[2][2] * tt[3][3] ) 
    + tt[1][2] * tt[1][3] * tt[2][3] + 1./27.; 
  double pTemp = max( min( rCoef / pow3(qCoefRt), 1.), -1.);
  double pCoef = cos( acos(pTemp) / 3.);
  double pCoefRt = sqrt( 3. * (1. - pow2(pCoef)) );
  eVal1 = 1./3. + qCoefRt * max( 2. * pCoef,  pCoefRt - pCoef);
  eVal3 = 1./3. + qCoefRt * min( 2. * pCoef, -pCoefRt - pCoef);
  eVal2 = 1. - eVal1 - eVal3;

  // Begin find first and last eigenvector.
  for (int iVal = 0; iVal < 2; ++iVal) {
    double eVal = (iVal == 0) ? eVal1 : eVal3;

    // If all particles are back-to-back then only first axis meaningful.
    if (iVal > 1 && eVal2 < EIGENVALUEMIN) {
      if (nBack < TIMESTOPRINT) 
      ++nBack;
      return false;
    }

    // Set up matrix to diagonalize.
    double dd[4][4];
    for (int j = 1; j < 4; ++j) {
      dd[j][j] = tt[j][j] - eVal;
      for (int k = j + 1; k < 4; ++k) {
        dd[j][k] = tt[j][k]; 
        dd[k][j] = tt[j][k]; 
      }
    }

    // Find largest = pivotal element in matrix.
    int jMax = 0;
    int kMax = 0;
    double ddMax = 0.;
    for (int j = 1; j < 4; ++j) 
    for (int k = 1; k < 4; ++k) 
    if (fabs(dd[j][k]) > ddMax) {
      jMax = j;
      kMax = k;
      ddMax = fabs(dd[j][k]);
    }

    // Subtract one row from the other two; find new largest element. 
    int jMax2 = 0;
    ddMax = 0.;
    for (int j = 1; j < 4; ++j) 
    if ( j != jMax) {
      double pivot = dd[j][kMax] / dd[jMax][kMax];
      for (int k = 1; k < 4; ++k) {
        dd[j][k] -= pivot * dd[jMax][k];
        if (fabs(dd[j][k]) > ddMax) {
          jMax2 = j;
          ddMax = fabs(dd[j][k]);
	}
      } 
    }

    // Construct eigenvector. Normalize to unit length. Random sign.
    int k1 = kMax + 1; if (k1 > 3) k1 -= 3;
    int k2 = kMax + 2; if (k2 > 3) k2 -= 3;
    double eVec[4];
    eVec[k1]   = -dd[jMax2][k2];    
    eVec[k2]   =  dd[jMax2][k1];    
    eVec[kMax] = (dd[jMax][k1] * dd[jMax2][k2]
      - dd[jMax][k2] * dd[jMax2][k1]) / dd[jMax][kMax];
    double length = sqrt( pow2(eVec[1]) + pow2(eVec[2])
      + pow2(eVec[3]) );
    if (gRandom->Uniform()> 0.5) length = -length;

    // Store eigenvectors.
    if (iVal == 0) eVec1 = TLorentzVector( eVec[1] / length,
      eVec[2] / length, eVec[3] / length, 0.);
    else eVec3 = TLorentzVector( eVec[1] / length,
      eVec[2] / length, eVec[3] / length, 0.);
  }

  // Middle eigenvector is orthogonal to the other two.
  eVec2 = cross3( eVec1, eVec3);
  //if (Rndm::flat() > 0.5) eVec2 = -eVec2;
if ( gRandom->Uniform() > 0.5) eVec2 = -eVec2;

  // Done.
  return true;

}

  // Return info on results of analysis.
  double sphericity()      const {return 1.5 * (eVal2 + eVal3);}
  double aplanarity()      const {return 1.5 * eVal3;}
  double eigenValue(int i) const {return (i < 2) ? eVal1 :
    ( (i < 3) ? eVal2 : eVal3 ) ;}
  TLorentzVector eventAxis(int i)    const {return (i < 2) ? eVec1 :
    ( (i < 3) ? eVec2 : eVec3 ) ;}

  // Tell how many events could not be analyzed.
  int nError() const {return nFew + nBack;}

private: 

  // Constants: could only be changed in the code itself.
    int    NSTUDYMIN, TIMESTOPRINT;
    double P2MIN, EIGENVALUEMIN;

// The scalar product of two three-vectors.

double dot3(const TLorentzVector& v1, const TLorentzVector& v2) {
  return v1.X()*v2.X() + v1.Y()*v2.Y() + v1.Z()*v2.Z();
} 

//*********

// The cross product of two three-vectors.

TLorentzVector cross3(const TLorentzVector& v1, const TLorentzVector& v2) { 
  TLorentzVector v; 
  v.SetX(v1.Y() * v2.Z() - v1.Z() * v2.Y());
  v.SetY(v1.Z() * v2.X() - v1.X() * v2.Z());
  v.SetZ(v1.X() * v2.Y() - v1.Y() * v2.X());
   v.SetT(0.);
return v; 
}

double max(double x1, double x2) const {if (x1>=x2) return x1;return x2;} 
double min(double x1, double x2) const {if (x1<=x2) return x1;return x2;}
double pow2(double x1) const { return pow(x1,2);}
double pow3(double x1) const { return pow(x1,3);}

  // Properties of analysis.
  double power;
  int     powerInt; 
  double powerMod;

  // Outcome of analysis.
  double eVal1, eVal2, eVal3; 
  TLorentzVector   eVec1, eVec2, eVec3; 

  // Error statistics;
  int    nFew, nBack;

public:

ClassDef(Sphericity,2)

};

//*********
 
// Constants: could be changed here if desired, but normally should not.
// These are of technical nature, as described for each.




#ifdef __MAKECINT__
ClassImp(Sphericity)
#endif



#endif